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A polymer chain tethered to a surface may be compact or extended, adsorbed or desorbed, 
depending on interactions with the surface and the surrounding solvent. This leads to a rich phase 
diagram with a variety of transitions. To investigate these transitions we have performed Monte 
Carlo simulations of a bond-fluctuation model with Wang-Landau and umbrella sampling algorithms 
in a two-dimensional state space. The simulations' density of states results have been evaluated for 
interaction parameters spanning the range from good to poor solvent conditions and from repulsive 
to strongly attractive surfaces. In this work, we describe the simulation method and present results 
for the overall phase behavior and for some of the transitions. For adsorption in good solvent, we 
compare with Metropolis Monte Carlo data for the same model and find good agreement between 
the results. For the collapse transition, which occurs when the solvent quality changes from good 
to poor, we consider two situations corresponding to three-dimensional (hard surface) and two- 
dimensional (very attractive surface) chain conformations, respectively. For the hard surface, we 
compare tethered chains with free chains and find very similar behavior for both types of chains. For 
the very attractive surface, we find the two-dimensional chain collapse to be a two-step transition 
with the same sequence of transitions that is observed for three-dimensional chains: a coil-globule 
transition that changes the overall chain size is followed by a local rearrangement of chain segments. 



I. INTRODUCTION 

Chain molecules anchored at surfaces are a part of 
many physical systems. Examples include polymer 
chains grafted to colloidal particles,^ proteins project- 
ing from cell membranes,'^ block copolymers at liquid air 
interfaceSjl^ and long chain molecules attached to inor- 
ganic surfaces for study with atomic force microscopes.'^ 
Parameters such as composition of the chains, grafting 
density, and interactions with the surface and the envi- 
ronment affect the properties of these systems by deter- 
mining the conformations of the chains. In this work we 
focus on the effects of solvent quality and surface interac- 
tions and investigate individual polymer chains tethered 
to a flat, impenetrable surface. A tethered chain near an 
attractive surface has many of the conformation charac- 
teristics of a free chain adsorbed from solution. Likewise, 
a tethered chain near a hard surface has much in common 
with an isolated free chain. These similarities allow us 
to compare with known results for free chains to validate 
our approach. They also suggest that our new results are 
relevant to the process of surface adsorption from dilute 
polymer solutions. 

As is well known, the collapse transition that a flexible 
polymer chain undergoes in dilute solution when the sol- 
vent quality deteriorates is one of the fundamental prob- 
lems in the statistical mechanics of polymers. It is of 
crucial importance in understanding the phase diagrams 
of polymer solutions and still offers surprising insights.^ 
Of similar importance is the adsorption transition in 
good solvent, where the conformations of a tethered chain 



change from "mushroom" to "pancake" configurations as 
the strength of the surface-monomer attraction increases. 
This problem has attracted longstanding intense atten- 
tion as a basic phenomenon of polymer chains interacting 
with interfaces.'^ Of course, both "transitions" are sharp 
thermodynamic phase transitions only in the (thermody- 
namic) limit when the chain length tends to infinity, but 
understanding how this limit is precisely approached is 
one of the challenges here. 

In the present work, we go one step beyond the prob- 
lems outlined above, by considering the interplay of ad- 
sorption and collapse. This problem is essential for un- 
derstanding polymeric aggregates on surfaces when no 
solvent is present( e.g. adsorbed polymers on a surface 
exposed to air!). It is also a first step towards the treat- 
ment of heteropolymers at surfaces, a system suited to 
give insight into the behavior of biopolymers at biologi- 
cal interfaces, membranes, etc. since these polymers of 
biological interest often assume rather dense conforma- 
tions. However, even the case of a homopolymer at a 
surface undergoing adsorption competing with collapse 
is a very difficult problem, and despite earlier numerical 
studies^ 9 10 ii 12 13 14 15 16 17 ^q^- j^^y understood. In 

particular, the question how the polymer model may af- 
fect the observed phenomena merits further investiga- 
tion. In most of the existing work, the polymer has 
been modeled as an interacting self-avoiding walk on a 
simple cubic lattice (ISAW). In the ISAW model, the 
bond length is fixed to the size of the lattice constant 
and the bond angles are restricted to the values of 7r/2 
and TT. This leads to very specific chain conformations in 
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poor-solvent conditions, where a chain near an attract ive 
surface fills the volume of a rectangular boxl ^^ l ^^ ^^ In 
contrast, an off-lattice model of a polymer chain under 
such conditions has chain conformations whose shapes 
resemble the spherical caps that are typical for partially 
wetting liquid dropsl^In this work, we employ the bond 
fluctuation (BF) model to investigate tethered chains. In 
the BF model, there are five allowed b ond leng ths leading 
to a large set of allowed bond angles i^^IMlI] terms of 
chain conformations, the BF model may be considered an 
intermediate between the ISAW lattice model and coarse- 
grained off-lattice models such as bead-spring models. 
In crystallizable homopolymers such as polyethylene and 
poly(ethylene oxide) the local structure of the chains 
leads to an alignment of chain sections and lamella for- 
mation upon crystallization.^^ The necessity for meso- 
scopic conformational ordering accompanying this align- 
ment gives rise to a large entropic barrier to crystal- 
lization and makes kinetic effects a dominant factor in 
polymer crystallization. Therefore it is difficult to sepa- 
rate thermodynamic effects from non-equilibrium effects. 
In coarse-grained models of polymers, such as the BF 
model and typical bead-spring models, one can choose 
force fields with sufficient local flexibility to allow for 
crystallization by local rearrangement of monomers. This 
reduces the entropic barrier to crystallization and allows 
the systems to reach ordered equilibrium states. While 
their lack of lamella formation hampers the compari- 
son with experimental data for simple crystallizable ho- 
mopolymers, the lower entropic barriers of coarse-grained 
models allow us to study equilibrium aspects of polymer 
crystallization. The BF model has proven to be a very 
useful model for a range of polymeric systems (see, for 
example, Refs. ^l9]20|2H23|24i25]) but has not yet been 
investigated in the present context. 

The combined effects of surface interactions and sol- 
vent quality lead to a variety of transitions for a poly- 
mer chain tethered to a surface. In Fig. [T| we present 
a schematic phase diagram that represents insights from 
theoretical and simulation work on lattice models of teth- 
ered polymer chains^ 9 lo ii 12 13 14 including results pre- 
sented in this work. The field variables /3s and /3f, are 
a combination of inverse temperature and interaction 
parameters^ and will be defined in Section III The 
field Ps describes the effect of surface interactions and 
increases with increasing attraction between the surface 
and the chain segments. The field /3t describes the net 
interactions between monomers of a polymer chain in sol- 
vent. Ph is small for a chain in a good solvent and in- 
creases with decreasing solvent quality, corresponding to 
increasing net attractive interactions between monomers. 

The regions marked DE (desorbed extended) and AE 
(adsorbed extended) correspond to good solvent condi- 
tions. With field values in the DE region, a tethered 
chain assumes extended three dimensional conformations 
that are often described as mushrooms. An increase of 
the surface field (3s leads to chain adsorption. For fields 
in the AE region, the chain conformations are (nearly) 
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FIG. 1: Schematic phase diagram for tethered chains of finite 
length in the space of field variables Pb and /3s (see Sec. Ill for 
the formal definition). As discussed in the text, increasing val- 
ues of /3s and /3i, correspond to increasingly attractive surface 
and monomer-monomer interactions, respectively. The lines 
indicate transitions between states identified by the following 
abbreviations: DE for desorbed extended (mushroom), AE 
for adsorbed extended (pancake), DC for desorbed compact, 
AC for adsorbed compact, and LS for layered states. The 
solid lines indicate transitions that are expected to become 
true phase transitions in the limit of infinite chain length. 
The dashed lines represent structural transitions observed for 
finite size chains only. The ending of the dashed lines indi- 
cates that, in simulations, these structural transitions can no 
longer be uniquely identified in regions of field parameters, 
where several transitions compete with each other. 



two dimensional and extended; they are sometimes called 
pancake conformations. The transition from DE to AE 
states has been investigated with theoretical and sim- 
ulation methods for a variety of models (see, for ex- 
ample, Refs. |18|26|27|28l29] ). In the limit of infinite 
chain length, the adsorption transition is sharp and cor- 
responds to a multicritical point. 

The regions marked DC (desorbed compact), AC (ad- 
sorbed compact), and LS (layered states) correspond to 
poor solvent conditions. Subject to fields in the DC re- 
gion, a chain assumes compact three-dimensional globule 
conformations and for fields in the AC region the cor- 
responding two-dimensional conformations. The region 
labeled LS is characterized by a competition between the 
effects of attractive surface and monomer interactions. 

For weak surface fields (3s, an increase in the field Pi, 
leads to chain collapse, i.e. the coil-globule transition 
from desorbed extended (DE) to desorbed compact (DC) 
states. For free chains, the coil-globule transition is a con- 
tinuous transition which may be followed by a first-order 
transition that is associated with spatial ordering of the 
chain segments (see, for example, Refs. |30I31I32] V From 
theoretical work on the ISAW model, the line describing 
the coil-globule transition is expected to be a straight 
line parallel to the Ps axi^ and perpendicular to the ad- 
sorption transition line. For strongly attractive surfaces. 
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i.e. for large /3s values, an increase in f3b leads to a transi- 
tion from adsorbed extended (AE) to adsorbed compact 
(AC) states. In analogy with the coil-globule transition 
of a free chain in two dimensions, this transition is ex- 
pected to be continuous for ISAW models.^ 

Adsorption in poor solvent occurs for high f3b values as 
the surface field /3s is increased. It takes a tethered chain 
from desorbed globule conformations in three dimen- 
sions (DC) to compact, single layer conformations (AC) 
through intermediate states (LS) that are not yet fully 
understood and appear to be model-dependent. Singh 
et a^l ^^ l ^^l^^ ] investigated untethered chains by exact enu- 
merations of interacting self-avoiding walks on a simple 
cubic lattice (ISAW model). In poor solvent conditions, a 
free chain assumes its most compact conformation which, 
in the ISAW model, has the shape of a cube with side 
length I (X iV^/^, where N is the length of the chain. In 
the presence of a (slightly) attractive surface, the cubic 
globule attaches itself to the surface and makes a num- 
ber of surface contacts proportional to the area of a face 
of the cube, i.e. cx N'^/^. Singh et aZ.ISIIIIlsl interpret 
the transition from the detached to the attached globule 
as a surface transition from the DC phase to a "surface 
attached globule" phase, which occupies the LS region 
in Fig. [T] and is followed by a second phase transition at 
high /3s values to the AC phase. For tethered chains in the 
ISAW model, theoretical arguments suggest that there is 
only one true phase transition (non-analyticity of the free 
energy in the limit of infinite chain length) as /3s is var- 
ied at constant /3h.'^More recent simulations of tethered 
chains in the ISAW model by Krawczyk et al^ found the 
chains to undergo a series of layering transformations in 
the LS region corresponding to transitions between com- 
pact chain conformations with different numbers of layers 
parallel to the surface. This work identified scaling re- 
lations which suggest that, with increasing chain length, 
the transitions become sharper and more numerous, and 
that the transitions merge into a single transition in the 
limit of infinite chain length. In the diagram of Fig. [T| 
we represent the true phase transition by a solid line be- 
tween the AC and LS regions. The dashed lines in Fig.[l] 
indicate structural/layering transitions that may not be- 
come phase transitions in the infinite chain limit. These 
transitions are important, however, since for many prob- 
lems the behavior of finite-length chains is of physical 
interest. 

In this work, we investigate transitions of tethered 
polymer chains with the aid of density-of-states simu- 
lations of a bond-fiuctuation model. In the work pre- 
sented here we focus on characteristic features of adsorp- 
tion in good solvent (DE to AE) and chain collapse in 
two and three dimensions (DE to DC and AE to AC). 
We also construct a phase portrait for our longest chain. 
In a forthcoming publication we shall discuss adsorption 
in poor solvent and investigate some states and transi- 
tions in detail. In this article, we introduce the bond- 
fiuctuation model for a tethered chain and describe the 
simulation methods in Sec. In] and the data evaluation 



in Section III Results of our simulations are presented 



and discussed in S ec.|IV| followed by a brief summary and 
conclusions in Sec.lVl Technical details of the simulations 
undertaken for this work are presented in Appendix [X] 



II. MODEL AND SIMULATION METHOD 

In this section, we present the model for a surface- 
attached polymer chain and the simulation methods em- 
ployed in this work. Details of the simulation protocol 
that are of interest to those who want to reproduce or 
extend the work but that are not required for an under- 
standing of the following sections have been relegated to 
Appendix 

In the bond fiuctuation (BF) model,'^^'^^!!!] monomers 
of a polymer chain occupy sites on a simple cubic lat- 
tice. The bond lengths between monomers are allowed 
to vary between b = 2a and b = vTOa, where a is the 
lattice constant, which we set to unity, a — I. A teth- 
ered chain is represented by a chain whose first monomer 
is fixed just above a hard surface. The position of the 
surface is the x-y plane at position z = in a Cartesian 
coordinate system; the coordinates of the fixed monomer 
are (1,1,1). Monomers at z = 1 are considered to be 
in contact with the surface and contribute an amount eg 
to the energy. The interactions between monomers de- 
pend on the distance. A pair of monomers i and j at a 
distance r^- contributes an amount of eb to the energy 
when 4 < r|j < 6. Distances rfj < 4 are prohibited by 
hard core repulsion, while monomers do not interact for 



r~j > 6. The total energy of the system is given by 



E{ns,nb) = Us^s + n-beb, 



(1) 



where Us and Ub are the number of monomer-surface and 
monomer-monomer contacts, respectively. We refer to a 
pair of contact numbers, (ris, nf,), as a state of the system. 

In the Monte Carlo simulations described here, two 
types of elementary moves were carried out; local moves 
that displace a monomer to a nearest neighbor lattice 
site and pivot moves about the z-axis. Simulations were 
performed for chains of length = 16, iV = 32, and 

= 64. In all simulations, N local move attempts were 
followed by ten pivot move attempts. We refer to this 
sequence as one Monte Carlo (MC) step. 



A. Wang-Landau algorithm for a 2-dimensional 
state space 

The density of states (dos), g{ns,nb) is the number 
of configurations for a given state (ris , n ^) of the sys- 
tem. The Wang-Landau (WL) algorithirPUMI jg an it- 
erative Monte Carlo simulation method for construct- 
ing the density of states. Originally formulated for a 
one-dimensional state space, characterized by the total 
energy, it has been extended to state spaces of higher 
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dimensions.l^^'^niSZI in this work, the state space is the 
two dimensional space of contact numbers and ni,. 

In the WL algorithm, an elementary move attempt 
from a state {ns,nb) to a state (n'^jTiJ^) is accepted with 
probability 



g{ns,nb) _ 



(2) 



where g(ns,n{,) is the current estimate for the density of 
states. 

The initial guess for the density of states is g{ns, Ub) = 
1 for all states {ns,nb). At each iteration level, elemen- 
tary move attempts are followed by an update of the 
density of states and the histogram of visits to the states, 
h{ns, Ub). After an attempted move from a state (n^, nb) 
to a state {n'^, the updates are given by 



if accepted: 
if rejected: 



ln(g«,nfc)) 
h{n',,n[) 

ln{g{ns,nb)) 
h{ns,nb) 



ln(5K,n',))+ln(/), 

ln{g{ns,nb)) + ln{f), 
h{ns,nb) + 1, 



(3) 



(4) 



where /, with / > 1, is the refinement factor. An it- 
eration is considered complete when the histogram sat- 
isfies a fl atness criterion. In a typical Wang-Landau 
simulation ,1^31211 a histogram is considered flat when the 
smallest value in the histogram is at least 80% of the 
average value of the histogram entries. At that point, 
the histogram is reset to for all states and the refine- 
ment factor is reduced before the next iteration is started. 
The adjustments to the density of states become smaller 
with each iteration level; we used refinement levels with 
ln(/fc) = 2-('=-i) for fee {2,20}. 

For tethered chains, the range of accessible states 
{ug, Ub) is not known a priori. Since the maximum num- 
ber of surface contacts ris is equal to the number of 
monomers N, the maximum value of the ratio Ug/N is 
unity for all chain lengths. However, the maximum value 
of Ub /N, the number of monomer-monomer contacts per 
monomer, increases with chain length as illustrated in 
Fig. |2] In addition, for each chain length, the maxi- 
mum number of monomer-monomer contacts decreases 
with increasing number of surface contacts. In the work 
presented here, we allow the range of considered states to 
change during the simulations. During the Wang-Landau 
simulations, when a configuration to a previously unvis- 
ited state (r7,|,,n'g) appears, the new state is assigned the 
initial dos value g{n'g,n'fj) = 1 and the move is accepted 
or rejected with the usual criterion. 

For simulations of free chains in the bond fluctuation 
model, it was found that updating after accepted moves 
only may lead to shorter simulation times without af- 
fecting the results for the density of states In the work 
presented here, we carried out simulations where either 
both updates, Eqs. ([3 1 and Q, or only the updates after 
accepted moves, Eq. (SI), were carried out. 
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FIG. 2: Ranges of accessible states for tethered chains of 
lengths TV = 7, 16, 32, and 64. The shaded areas approximate 
the ranges, the symbols represent the realized states with the 
highest number of bead contacts nt for given number of sur- 
face contacts n^. The results for N — 7 were obtained in 
exact enumeration. 



B. Multiple replica algorithm 

In a multiple replica simulation, the state space is 
divided into a number of overlapping regions. Stan- 
dard simulation steps are carried out separately in each 
replica. After each MC step, an attempt is made to ex- 
change the chains of two neighboring replicas if their 
configurations belong to states in the overlap region of 
the replicas. The acceptance criterion for such a replica- 
swapping move is 



p{[{ns,nb)i,{n'g,n'^)^+i] {n'g,n'^)i, {ns,nb)^+l]) 

g{ns,nb)i+i (5) 
g{ns,nb)^ ' 



where i and i + 1 are the indices for neighboring replicas. 
When using the method to generate the density of states, 
the dos values are updated after each attempted move in 
the usual way, see Eqs. ([3| and Q. The density of states 
for the whole state space is obtained at the end of the 
simulation by combining the results for each replica. For 
the results presented here, we used regions that were de- 
fined by intervals in the number of bead-bead contacts Ub 
and contained all possible surface contacts. We experi- 
mented with regions of different size and different ranges 
of overlap and found that small but systematic errors in 
the density of states occurred unless each state was cov- 
ered by the same number of replicas. We worked with a 
total of three to five overlapping replicas covering each 
region with two replicas. 
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C. Global update algorithm 

When the WL algorithm is appHed to large two- 
dimensional state spaces much time is spent accumulat- 
ing the large ln{g) values belonging to the interior of the 
state space. Introducing global update^^ improves the 
efficiency by allowing the simulation to spend more time 
exploring the edges of the current state space. In a global 
update algorithm, a simulation is started at a given re- 
finement level and the current values of ln(g(rts, rib)) are 
compared with a threshold value uj. Once density of 
states values above the threshold, lii(g(ns,nh)) > lo are 
found, they are augmented in the following way: 

ln(5) ^ ln(g) + ^cxp (^^^^) ^iHo) " (6) 

where is the Heaviside step function. The exponen- 
tial function dampens the shift from a maximum value 
of K for states with ln{g) values well above the thresh- 
old to zero for states at the threshold; the parameter A 
determines the range over which the shift is phased out. 
The simulation then continues with standard steps and 
local updates to the density of states. Immediately after 
a global shift, states that were outside the shifted region 
are preferentially sampled. With time, their hi{g) val- 
ues increase thereby decreasing the acceptance rates for 
moves to those states and leading to uniform sampling 
over previously shifted and unshifted states. Once uni- 
form growth over the current state space has resumed the 
next global update is considered. This continues until a 
flatness criterion is satisfied or the simulation is termi- 
nated by hand. 

To determine when uniform growth has been achieved 
in our simulations, we monitor the difference Aln(g) be- 
tween the current value of the log-density of states and 
that at the time of the global shift. When for each previ- 
ously shifted state the difference A ln{g) is larger than a 
fixed parameter p, we conclude that uniform growth has 
been achieved. We found that the values for the param- 
eters K, A, u and p can have a large effect on the quality 
of the results. In particular, systematic deviations in the 
results for the density of states occur when the global 
updates are too large and occur too frequently. 



D. Umbrella sampling 

In multicanonical or umbrella sampling simulations 
with single histogram reweighing, a good estimate for the 
density of states g{ns,nh) is required as input (see, for 
example, Refs. [21|39)40|41j ). We performed umbrella 
sampling simulations with a WL result as the initial dos 
for all simulation results reported here. During the sim- 
ulation, elementary moves are accepted with the usual 
criterion, see Eq. ([2]), and the histogram is updated after 
accepted and rejected moves. The dos values, however, 
are updated only at the end of the simulations. In a 



typical single histogram reweighing step, the logarithm 
of the final histogram entry h(ns,nb,tf), is added to the 
original log-density of states values, 

ln(.g'(7is,7ib)) = \n{g{ns,nb)) + \n{h{ns,nb,tf)). (7) 

In the work presented here, we record cumulative 
histograms at regular intervals, determine the slope 
m{ns, Hh) of the histogram entries as a function of time, 
and estimate the final histogram entry from the slope 
h{ns,rn,,tf) ~ m{ns,'nb)tf . This gives a slightly more 
reliable dos estimate for states that are not visited in 
every block. 

E. Metropolis algorithm 

For comparison with results from the density of states 
algorithms and in order to sample efficiently some parts 
of phase space during the production stage, some sim- 
ulations were performed with a Metropolis acceptance 
criterion at fixed fields (3^ and (3},. The probability for 
accepting a move from a state (n^, rib) to a state (n^, 
in the Metropolis algorithm may be written as 

p((n,,nb) ^ 

= min {e^])[fis{n'^ - n^) + /3b{n'b - rib)], 1) • 

F. Production stage 

In a production simulation, chain conformations are 
evaluated at regular intervals to accumulate configura- 
tional properties as a function of the pairs of contact 
numbers {ns,nb). We determined a range of configu- 
rational properties including chain dimensions and den- 
sity profiles. Furthermore, configurations were stored 
for a detailed analysis of the chain structure. In this 
work, we present results for the squared bond length 
= ^■(ri_|_i—ri)^/(A^—l), where is the position vec- 
tor of monomer i, and for the radius of gyration Rg. In 
order to investigate the effect of the surface, we calculate 
parallel and perpendicular contributions to Rg according 
to 

i<j 

where i,j € {!,..., N} and where Xi,yi, are the Carte- 
sian coordinates of monomer i, and R^, = Rj, , + R?„ 

When production simulations are performed with the 
acceptance criterion of the Wang-Landau algorithm all 
states (rig, rib) are visited approximately equal numbers 
of times. This may lead to insufficient sampling of states 
(ris, rib) with a diverse set of chain conformations. In or- 
der to improve the statistics for the configurational prop- 
erties belonging to such states (n^, rib), we performed ad- 
ditional production simulations with a Metropolis Monte 
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Carlo algorithm for fields between /?s = /?(, = and 
j3s = I3b = —2. In these simulations, the visitation his- 
tograms are approximately Gaussian with maxima at the 
most probable states for these fields. 

For all chain lengths, we performed more than one 
type of produ ction simulation. During production, block 
averagea^IEll were calculated and the statistical uncer- 
tainty of an average value was determined from the stan- 
dard deviation of the block averages. Results from dif- 
ferent production simulations were combined as weighted 
averages, where the weights were based on the uncer- 
tainty estimates obtained for each quantity in the indi- 
vidual production simulations. 



III. EVALUATION OF THE DENSITY OF 
STATES 

The canonical partition function for a tethered chain 
is given by 



(11) 



where (3 — T is the temperature, fee is Boltz- 

mann's constant, and g{ns,ni,) is the density of states. 
The energy of the state, E{ns,nb) = UgCs + nbtb, de- 
pends on the contact numbers Ug and Ub and the inter- 
action parameters and tb- In this work, we employ 
combinations of the interaction p aram eters and the in- 
verse temperature as field variabled^^ll 

A - Pb = (12) 

and write the partition function as 

Z(/3„/3b) = ^ 5(n„n,)e^="=e^'"^ (13) 

For fixed fields (3^ and (3b, the probability P{ns, rib, (3s, /?&) 
for states with contact numbers (n^, Ub) is given by 

P{ns,nb;Ps,(3b) = ^g{ns,nb)el'^"^ e^'^^ (14) 

and the average values of quantities Q{ns,nb) at given 
fields (3s and (3b are calculated from 

Q{ns,nb)P{ns,nb;(3s,(3b) = Q{Ps,(3b)- (15) 

At the end of this section, we describe how we estimate 
uncertainties for the averages determined in this way. 
With the definition of a free energy 



Gi(3s,Pb)^-HZ), 



(16) 



the average number of surface contacts (rig) and the av- 
erage number of monomer contacts {ub) are conjugate to 



the field variables Ps and Pb, respectively, 

dG \ 
dPs),,' 



{nb) 



(17) 
(18) 



The second derivatives of the free energy with respect to 
the fields define susceptibilities Xfu^i M; ^ ^ {-^j ^} that are 
related to the fiuctuations in the number of surface and 
monomer contacts 



Xhb 



Xsb 



<fiG_ 



d^G 

dPsPb 



{Ub 



[nsUb) 



>{nb) 



(19) 
(20) 
(21) 



In terms of these susceptibilities, the heat capacity of the 
system is given by 

C = P^ {(E^) - (E)^) = P^,Xss + Phbb + 2PsPbXsb. (22) 

There are alternative methods to evaluate the density 
of states which we would like to discuss briefly. The den- 
sity of states g{ns,nb) provides the entropy as a function 
of the contact values up to an arbitrary constant 5*0 



S{ns,nb) fee \n{g{ns,nb)) + Sq. 



(23) 



Starting from the entropy as the thermodynamic poten- 
tial, the fields and susceptibilities may be expressed in 
terms of derivatives of the log density of states. For 
example, the fields may be obtained as Ps{ns,nb) = 
-{dln{g)/dns)nt and Pb{ns,nb) = -{dln{g)/dnb)n,- In 
the thermodynamic limit, N oo, results from differ- 
ent statistical ensembles are equivalent, however, for our 
chain lengths finite size effects may play a role. Since an 
evaluation of the entropy is hampered by the need to take 
numerical derivatives of discrete variables, we employ the 
free energy G{Ps,Pb) unless otherwise indicated. 

In some cases, it is convenient to work with a ther- 
modynamic potential that is a function of the number 
of surface contacts Us and the bead-contact field Pb- 
This free energy A(ns,Pb) may be considered a Legendre 
transforrrP^ of either the entropy S or the free energy G 
and may be written as A(ns,Pb) = ^Iii(Za) , with 



Z^(n„/3fc) = ^5(n„nb)e'3^"'' 



(24) 



In this case, the probability for a state with nb bead 
contacts is given by 



1 



Pinb;ns,Pb) = ^g{ns,nb)e'''''\ (25) 
Za 



7 



which allows us to calculate quantities for given and 
f3b without taking numerical derivatives. For example, 
the average number of bead contacts is determined from 



(3b). 



(26) 



In general, there are two sources of error when average 
quantities (Q) — Q{(3s,(3b) for given field values f^s and 
(3b are obtained from the contact number dependent val- 



ues Q(ris, Ub) with the aid of Eq. ( 15 1. The first source of 
error is the uncertainty in the density of states, this is the 
only uncertainty for quantities such as the susceptibili- 
ties, see Eqs. (19 1 - (21 1. To estimate this uncertainty. 



we perform several simulations for the same chain length 
and calculate the mean value lii(g(ns, nb)) and the stan- 
dard deviation <J\n(g){ns,nb) for each state {ns,nb). The 
second source of error is the uncertainty in production 
results, which we obtain from block averages. For prop- 
erties such as the chain dimensions, see Eqs. 
this is typically the larger source of error. 

In order to estimate the effect of the uncertainties in 
the density of states on calculated properties Q, we gen- 
erate five synthetic densities of states \^{gi) by draw- 
ing randomly from Gaussian distributions centered on 
\n{g{ns,nb)) with standard deviation crin(g) (n^ , rif,) for 
each state {ns,nb). Average values {Q)i are calcu- 
lated from each synthetic density of states i and the 
mean Q and standard deviations crg.g of the results 
are determined. The uncertainties of production data 
are propagated through Eq. (15), evaluated with ln(g), 
to obtain the second contribution to the uncertainty 
CTQ^p. The combined uncertainty is estimated from ctq — 



Q,g ^ "Q,p- 

The uncertainty estimates, crg.g for results that derive 
their errors only from the density of states and ctq for 
results that involve production data, are shown as error 
bars in the figures of this work. Error bars are omitted 
when they are smaller than the symbol size or the line 
thickness. 



IV. RESULTS AND DISCUSSION 

A. Density of states 

Density-of-states results were obtained for chains of 
length iV = 16, = 32, and A^ = 64 with the simula- 
tions described in Section |IT] Figure [2] shows the ranges of 
accessible states for these chains and also includes results 
from an exact enumeration for chains of length N = 7. 
In each case, the maximum number of bead-bead con- 
tacts decreases with increasing number of surface con- 
tacts since chain conformations cannot maximize surface 
and bead contacts simultaneously. This leads to a com- 
petition between bead and surface contacts at high field 
values which increases with increasing chain length. 




FIG. 3: Density of states for tethered chains of length A = 32 
(top) and N = 64: (bottom). The surfaces (small symbols 
connected by straight line segments) represent the log density 
of states values, \n{g), as a function of surface contacts, Ua, 
and bead-bead contacts, nt- 



In Fig. [3] we present results for the log density of states 
for chains of length A^ = 32 and A^ = 64 (the results for 
A = 16 show the same qualitative behavior); some char- 
acteristic numerical values are presented in Table |T] For 
each surface contact value, Ug, and for small contact val- 
ues Hb, the density of states increases with rtf,. As nb 
increases further, ln(g) passes through a maximum and 
then decreases with increasing Ub- As discussed at the 
end of Sec. Ill the slope of \n{g{ns, rib)) at fixed Ug is re- 



lated to the expectation value of the field (3b{ns,nb). This 
suggests that states with small bead contact values are 
predominantly populated for negative (3b fields (repulsive 
interactions between the beads) while states with large 
bead contact values are populated for positive (3b fields 
(attractive bead-bead interactions). This is confirmed 
by an evaluation of the average bead contact number ac- 
cording to Eq. (26 1. In Fig. |4]we present results for the 
average number of bead contacts (nf,(ns, /5b)) as function 
of the number of surface contacts for five fields (3b- For 
each Us and f3b we also evaluated cumulative probabili- 
ties X^n^lo P{''T'si (3b', nb) and determined the bead contact 
numbers n*b,i and n*b,2 where the cumulative probabil- 
ity first exceeds 1/6 and 5/6, respectively. The states 
between n*f,,i and n*b.2 have significant probability of oc- 
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FIG. 4: Most probable states for given number of surface con- 
tacts Us and field Pb- The lines represent the average number 
of bead contacts {nb{ns, Pb)) as function of the number of 
surface contacts for five fields /3b, as indicated in the figure. 
The shaded areas surrounding the lines indicate the states 
that have significant probability of occupation; the sum of the 
probabilities associated with these states is 2/3. The symbols 
at the upper boundary of the graph indicate the maximum 
values of nj, for given n^. 



cupation and are indicated by the shaded areas in Fig. |4] 
As expected, the range of significantly populated states 
shifts from low-nt, to high-n;, values with increasing field 
Pb- For the largest field shown, f3b = 2, states on the 
upper rim of the range of accessible states start to be- 
come populated. Since the density of states values for 
rim states have a relatively large uncertainty, we restrict 
ourselves in this work to fields /3fc < 2. The width of 
the shaded areas in Fig. |4] is an indication for the size 
of the fluctuations about the mean value. These fluctua- 
tions are larger for /3b = 1 than for the other flelds, since 
/3b = 1 is closer to the fleld value of the coil-globule tran- 
sition (see Sec. IV C I than the others. Geometrically, the 



size of the fluctuations is associated with the inverse of 
the curvature of the log-density of states surface. 



For most accessible states, the density of states is a 
monotonously decreasing function of at constant ni,, 
i.e. the states are predominantly occupied for attractive 
surface interactions. Only for a range of states with small 
Us and very large Uf, values is the slope of ln(g) for fixed 
rib positive. These states become populated predomi- 
nantly when Ps is negative and /3b is large and positive. 
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FIG. 5: Average number of surface contacts per monomer 
(ris) /N and surface contact fluctuations, Xs/N, as a function 
of the surface field Ps for good solvent conditions {pb = 0). 
The lines represent results from the evaluation of the density 
of states for chains of length A*' = 64 (solid), 32 (dashed), and 
16 (dash-dotted), respectively. The graphs for {na)/N are 
monotonously increasing, those for Xs/N have a maximum in 
the transition region. 



B. Adsorption in good solvent 

To explore adsorption in good solvent conditions, the 
bead contact field is set to /3b = so that chain seg- 
ments interact with each other only through excluded 
volume interactions. Density of states and production 
results are evaluated as described in Sec. |III| for a range 
of surface fields /3s. In Fig. [5] we present results for the 
average number of surface contacts and their fiuctuations 
as a function of the surface field Ps- For negative Ps val- 
ues, the surface is repulsive and the number of surface 
contacts is near 1, the smallest possible value. As Ps in- 
creases, the surface becomes increasingly more attractive 
and the number of surface contacts increases until it ap- 
proaches N, the largest possible value. The transition 
between desorbed and adsorbed states becomes steeper 
with increasing chain length which is an expected finite 
size effect. The fluctuations Xs sn'e small for repulsive and 
very attractive surfaces and have a maximum in the tran- 
sition region. The location of the maxima are Ps = 1.29, 
1.24, and 1.19 for chains of length N = 16, 32, and 64, 
respectively. We use maxima in Xs to identify the tran- 
sition fields for surface adsorption in Sec. |IVE[ 

The shape of a tethered chain may be described with 



the contributions Ri 



to 



and Rl^y ^ 2Rl^^ 

the radius of gyration i?^ .12^ For repulsive or weakly ad- 
sorbing surfaces, the chain assumes mushroom configu- 
rations whose extensions parallel and perpendicular to 



the wall are comparable, Rg.i^ 
becomes more attractive, Rg.±_ 



> R 



v,\\- 



As the surface 



decreases and R„ 



creases, so that the ratio Rg z/^g 



0-5Rl^ 



IK 



m- 
de- 



creases rapidly in the transition region. A scaling analy- 
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FIG. 6: Ratio R^.z/R^.xy of perpendicular and parallel con- 
tributions to the square radius of gyration as a function of the 
surface field /3a for good solvent conditions (/3t. = 0). The lines 
represent results from the evaluation of the density of states 
and production data for chains of length A'' = 64 (solid), 32 
(dashed), and 16 (dash-dotted); the inset shows an enlarge- 
ment of the region where the lines cross. For clarity, error 
bars for our calculated values are shown only in the inset. 
(The error bars generally increase with increasing A'^ and de- 
creasing /3s; for /3s = they are about twice as large as for the 
/3s-range of the inset). The filled symbols with error bars rep- 
resent Metropolis Monte Carlo results from this work. The 
open symbols represent Metropolis Monte Carlo results for 
chains of length = 20, 40, and 80 by Descas et al.W^ 



sis shows that the ratio is independent of chain length at 
the adsorption transition. In Fig. |6]we present results 
for the ratio Rl^z/^l,xy a function of (3s for /3b — 0. 
The lines represent results from an evaluation of our den- 
sity of states and production results for chains of length 
N — 16, 32, and 64. For comparison, we have also in- 
cluded Metropolis Monte Carlo results obtained by us for 
chains of length = 32 and 64, and by Descas et al:'^^ 
for chains of length N — 20, 40, and 80. The agreement 
between data obtained with different methods is good. 
While longer chains are required for a detailed compar- 
ison with scaling predictions, our results for z/-Rg xy 
show reasonable behavior. The lines for A^ = 64, 32, and 
16 intersect each other at about /3s = 0.96, as shown in 
the inset. Our value for the intercept, (3s = 0.96 ±0.02, is 
consistent with the results of Descas et al.,^^ who deter- 
mined an intercept of (3s — 0.98 ± 0.03 from data which 
included chains up to length 200. 



C. Chain collapse for a hard surface 

A polymer chain undergoes a transition from extended 
to compact conformations, the coil-globule transition, 
when the solvent quality changes from good to poor. For 



O 




FIG. 7: Heat capacities per monomer, C/N, as a function 
of the reduced temperature T* — 1/(3^. The solid (A'' = 64), 
short-dashed (A' = 32), and dash-dotted (A'' = 16) lines repre- 
sent results for chains tethered to a hard surface {(3s = 0). The 
long-dashed (A'^ = 64) and dotted (A'' = 32) lines represent 
results for free chains.^-^-^The estimated uncertainties of the 
results for the tethered chains are smaller than the line thick- 
ness, except for A'^ = 64 at very low temperatures, T* < 0.7, 
where they correspond to about twice the line thickness. 



the bond-fluctuation model with attractive bead-bead in- 
teractions, Ef, < 0, this transition may be induced by 
reducing the temperature. The continuous coil-globule 
transition is followed by a discontinuous crystallization 
transition at lower temperature. For the model we are 
studying, which has a range of attractive interaction 
< \/6, the difference between the transition temperatures 
decreases with chain length and, in the limit N — > oo, the 
transitions merge into a single first-order transition.l^Sl 

To explore chain collapse for tethered chains, we con- 
sider first a hard surface and set (3s = 0. In Fig. [7] we 
present results for the heat capacity, C, as a function of 
the reduced temperature T* — l.Q/(3h obtained by eval- 
uating our density of states results for chains of length 
N = 16, 32, and 64. For comparison, we also present re- 
sults for the heat capacity of free chains of length A^ = 32 
and A^ = 64, ob tained from Wang-Landau simulations 
of the BF model.l^^lMI por each chain length considered 
here, the heat capacity shows a peak at low tempera- 
tures followed by a shoulder at higher temperatures. For 
free chains, the peak is associated with the crystallization 
transition and the shoulder marks the coil-globule transi- 
tion; we find the same to be true for the tethered chains 
In Fig. 11 (a) we present an example for a compact and 
highly ordered conformation of a chain of length A^ = 64 



that is typical for the lowest temperatures. 

While the heat capacity values for the tethered and 
free chains coincide at high temperatures, there are some 
differences at lower temperatures. For chains of length 
A^ = 64 and temperatures smaller than T* w 1, the free 
chain heat capacity values are lower than those of the 
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FIG. 8: Monomer-monomer contact fluctuations, Xb/N, as a 
function of the field Pt for a hard surface, /3a = 0. The lines 
represent results from the evaluation of the density of states 
for tethered chains of length A*" = 64 (solid), 32 (dashed), and 
16 (dash-dotted). 
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FIG. 9: Monomer-monomer contact fluctuations, Xh/N, as a 
function of the field fib for a very attractive surface, /3s = 8. 
The lines represent results from the evaluation of the density 
of states for tethered chains of length = 64 (solid), 32 
(dashed), and 16 (dash-dotted). 



tethered chain. This is due to a difference in the sim- 
ulations; for the free chain, a cutoff of 272 bead con- 
tacts was used while the tethered chain results include 
bead contacts up to nh^max = 288. For somewhat higher 
temperatures, I T* < 1.5, there is also a small tem- 
perature shift (barely visible in the graph) between the 
heat capacity curves of the free and tethered N — 6A 
chains. For chains of length N = 32 and temperatures 
smaller than about T* = 1.2, i.e. below the coil globule 
transition, the heat capacity curve for the free chain is 
shifted to lower temperatures compared to that of the 
tethered chain. A comparison of the normalized den- 
sity of states results for the free and tethered chains of 
length iV = 32 shows that the number of available con- 
formations decreases with increasing bead contact num- 
ber more rapidly for free chains than for tethered chains. 
This suggests that the hard surface reduces the relative 
number of extended conformations more than the relative 
number of compact conformations. Since the fraction of 
monomers belonging to the surface of a conformation de- 
creases with chain length, it seems reasonable that the 
effect is larger for shorter chains; in our simulations, the 
effect has all but disappeared for N — 64. In the limit 
N — > oo, we expect the transition temperature for the 
three-dimensional coil-globule transition to be the same 
for free and tethered chains. 

To determine the location of the transitions it is con- 
venient to consider the fluctuations in the number of 
bead contacts Xb- These are related to the heat capacity 
through Eq. (22). In Fig. [s] we present tethered chain re- 
sults for Xb as function of the field Pb for chains of length 
N = 16, 32, and 64 at = 0. The graphs for iV = 32 
and N = 64: show two well separated maxima while the 
graph for iV = 16 has one maximum associated with the 
coil-globule transition and a shoulder associated with the 



globule-globule transition. The locations of the maxima 
for the coil-globule transitions are 0.91, 0.79, and 0.71 
for chains of length N = 16, 32, and 64, respectively. 
The field values for the globule-globule transitions are 
1.47 (estimated), 1.05, and 0.850 for N = 16, 32, and 
64 respectively. The values for the globule-globule tran- 
sitions compare well with those obtained from the peaks 
in the heat capacity, which yield (3b values of 1.48, 1.08, 
and 0.853 for N = 16, 32, and 64, respectively. We use 
maxima in Xb to identify the transition fields for chain 



collapse in Sec. IV E 



The transitions of single untethered polymer chains 
have been investigated intensivel y with a variet y of off- 
lattice polymer models.^" ^4 4 5|4B|47|48 | 49|5U|51 ] Except 
for the shortest chains and the most short-ranged inter- 
actions, the heat capacity as a function of temperature 
shows two prominent features; a peak at higher temper- 
atures that indicates the transition from disordered coil 
conformations to disordered globule conformations, and 
a peak at low temperatures that indicates the transition 
from a liquid-like disordered globule to a solid-like or- 
dered globule. Typically, the peak associated with the 
order-disorder transition, which is sometimes referred to 
as the "freezing" transition, sharpens with increasing 
chain length and is accompanied by a bimodal proba- 
bility distribution as expected for a discontinuous transi- 
tion. Our results for tethered and free chains show that 
the chain collapse in the bond fluctuation model is anal- 
ogous to chain collapse in off-lattice models. 



D. Chain collapse in two dimensions 

If a chain is completely adsorbed to a surface, a col- 
lapse transition between two-dimensional states may be 
observed. To investigate this transition, we chose a 



11 




0.0 0.5 



1.0 



1.5 



2.0 



6 







N = 64 


|dRVdPbl/N 


Ps-8 - 

✓ \ 

r \ 

y '\ 


-- IdB^/dPbl 
Xb/N 




A. 




\\ ^ 
~ \ ^ 











0.0 0.5 



1.0 1.5 

Pb 



2.0 



FIG. 10: Chain collapse for large surface field. The top panel 
shows the average square bond lengths B^, radius of gyration 
fig divided by the chain length A'^, and the average number of 
monomer-monomer contacts per bead (nt,) jN as a function 
of the monomer-contact field /S;, for a chain of length A?^ = 64 
and a surface field fia = 8.0. The bottom panel shows how 
these quantities change as the chain undergoes the collapse 
transition. The dashed and dash-dotted lines represent ab- 
solute values of the numerical derivatives {dR^g/d(3b)i3s ^nd 
{dB'^ /dl3t)i3^, respectively. The solid line represents Xb/N as 
in Fig. [I 



strong surface field, /S^ = 8, and evaluated the density of 



states as usual, with probabilities given by Eq. (14 1. For 
comparison, we also performed calculations using only 
density of states results with the maximum number of 
surface contacts, — N , and evaluated them with prob- 
abilities given by Eq. (251. For the range of f3b values 



relevant to the collapse transition, the results are indis- 
tinguishable. (For much higher values, the transition 
to layered states appears in the evaluation of the full 
density of states but not the restricted set.) In Fig. [9] 
we present results for the bead contact fluctuations Xb of 
tethered chains as a function of (3b for chains of length 
N = 16, 32, and 64 at (3^ = 8. For N = 16, only a single 
maximum is discernible, while for = 32 the maximum 
is broad and accompanied by a shoulder on the high /3f, 
side. For N = 64, a narrower peak in Xb is preceded by 
a shoulder on the low (3b side. In analogy with the two- 
stage collapse transition for chains tethered to a hard 




FIG. 11: Examples for compact conformations of chains of 
length A'^ = 64. (a) A highly ordered three-dimensional con- 
formation representative of the desorbed compact (DC) re- 
gion of the phase portrait in Fig. 12 (b) A highly ordered 
two-dimensional (single layer) conformation representative of 
the adsorbed compact (AC) region of the phase portrait. In 
both diagrams, the size of the circles corresponds to the hard 
core diameter of the beads; the bonds are shown as wide lines. 



surface, we expect the features at low and high (3b values 
to indicate coil-globule and globule-globule transitions in 
two dimensions, respectively. 

To investigate this further, we calculated the square 
radius of gyration and bond lengths of the chains. 
In Fig. [To] we present results for these quantities and 
their derivatives with respect to /3f, as a function of (3b 
for a chain of length = 64 in the surface field (3s — 8- 
For comparison, the average number of bead contacts, 
{nb)/N, and their fluctuations, Xb/N, are also shown. 
The graphs show that the largest changes in the radius of 
gyration occur near the shoulder of the Xb graph, while 
the largest changes in the bond length occur near the 
maximum of Xb- Since the radius of gyration measures 
the overall size of a chain conformation, while the bond 
lengths represent a small-scale property, the results sug- 
gest that a coil-globule transition is followed by a local 
rearrangement of segments. As in the three-dimensional 
case, the states beyond this latter transition are highly 
ordered, i.e., we interpret this "globule-globule transi- 
tion" as a transition from a fluid-like to a crystal-like 
state. In Fig. 11 (b) we show an example for such a 
highly ordered two-dimensional conformation of a chain 
of length N = 64. As expected for an order-disorder 
transition, we find a bimodal probability distribution for 
the occupation of states at the globule-globule transition. 

A two-stage transition from a disordered extended coil 
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FIG. 12: Phase portrait for a tethered chain of length A'^ = 64 
in the space of the surface field (/3s) and bead-contact field 
{f3b)- The regions are named as in Fig. [l] The solid lines 
represent maxima of surface (xs) and bead-contact (xt.) fluc- 
tuations, as explained in the text. The dashed lines are an 
estimate for the location of the coil-globule transition from 
the "shoulder" on the susceptibility Xb- The dotted lines rep- 
resent shallow maxima in the susceptibility Xa that depend 
sensitively on the details of the available compact chain con- 
formations. In the shaded area near the center of the diagram, 
the susceptibility "landscapes" are too complex to identify all 
of the maxima clearly. This is why some of the lines end 
rather than merge with other lines. 



through a disordered globule to an ordered globule in two 
dimensions has also been observed in simulations with 
a parallel tempering algorithm of an off-lattice (bead- 
spring) model.'^ The symmetry of the two-dimensional 
ordered conformations is model dependent. For the BF 
model employed in this work, the unit cell for the or- 
dered interior of Fig. 11 (b) has basis vectors ai — 2ay 
and a2 = ay -f 2ax, where x and y are the unit vectors 



in the x and y direction, respectively, and a is the lat- 
tice constant. The number of nearest neighbors (six) is 
the same as for a hexagonal lattice, which is the sym- 
metry for the two-dimensional ordered state of the off- 
lattice model.l^For the ISAW model, on the other hand, 
compact two-dimensional chain conformations have the 
symmetry of the square lattice with four nearest neigh- 
bors. It is an open question to what extent these kinds of 
two-dimensional states carry over to chemically realistic 
models including rough substrates. 



E. Phase portraits 

In the preceding sections we have employed maxima in 
the susceptibilities Xs and Xb to identify transition fields 
f3s and Pt for chain adsorption in good solvent and for 
chain collapse near hard and very attractive surfaces. In 
this section, we extend this approach to a wide range 
of conditions and construct a phase portrait for chains 



of length N = 64, shown in Fig. [12] that may be com- 
pared with the phase diagram of Fig. [l] discussed in the 
Introduction. 

From the maxima and "shoulders" of the susceptibil- 
ities we identified transition lines in the space of field 
variables /3s and (3b that separate the desorbed extended 
(DE), adsorbed extended (AE), desorbed compact (DC), 
adsorbed compact (AC), and layered state (LS) regions. 
The shaded area near the center of the phase portrait 
in Fig. [12] marks a range of field values where the sus- 
ceptibility "landscapes" are too complex to identify all 
of the maxima clearly. This is the reason why some of 
the lines end rather than merge with other lines. In this 
discussion, we focus on field values outside the shaded 
area. 

In the diagram of Fig. 



12 the horizontal solid lines 



separating the DE and DC regions represent maxima of 
Xb associated with the dual collapse transition in three 
dimensions. Similarly, the horizontal solid and dashed 
lines separating the AE and AC regions represent max- 
ima and "shoulders" of Xb associated with the collapse 
transition in two dimensions. For both collapse transi- 
tions, the fiuctuations in the number of surface contacts, 
Xs, are not significant. The vertical solid line separating 
the DE and AE regions represents maxima in Xs asso- 
ciated with adsorption in good solvent. For this transi- 
tion, the fiuctuations in the number of bead contacts, Xbi 
are small. The two solid lines that separate the regions 
of adsorbed states (AE and AC) from fully or partially 
desorbed states in poor solvent conditions (DC and LS) 
represent maxima of both surface and bead-contact fiuc- 
tuations. The dotted lines in the LS region represent 
shallow maxima of Xs- The error bars show the spread 
of the location of these maxima when different sets of 
density of states values are evaluated. The transitions 
associated with these maxima are structural, discontinu- 
ous transitions that depend sensitively on the details of 
the available compact chain conformations. 

The AC states for high Ps and f3b values have strictly 
two-dimensional, single-layer conformations with almost 
hexagonal symmetry as shown in Fig. 11 (b). The chain 
conformations for slightly lower (3s values, in the region 
between the two solid lines near the upper right corner of 
Fig. [12] have exactly two layers, each with nearly hexag- 
onal symmetry and shifted against each other by one lat- 
tice constant a of the underlying cubic lattice. The tran- 
sition from single to double layers is similar in nature to 
the layering transition observed in the ISAW model. In 
the limit of infinite chain length, the slope of the line 
separating the two-layer from the one-layer conforma- 
tions may be determined by equating the energy per seg- 
ment of the single and double layer conformations .1^^ For 
the ISAW model, the value of this slope is 1/2. For the 
bond-fiuctuation model employed in this work, the cal- 
culated value for infinitely long chains is 1/3 while finite 
size effects increase it to a value of 0.38 for the chains of 
length N = 64. States in the region marked LS in Fig. [12] 
are very different from the single and double-layer states 
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discussed above. The chain conformations are three- 
dimensional with a symmetry group corresponding to 
a cubic lattice with basis. The number of surface con- 
tacts decreases stepwise through the transitions marked 
by dotted lines until only a single surface contact is left 
in the region marked DC. 

For the model of a tethered chain investigated in this 
work, all transitions from good to poor solvent conditions 
are dual in nature. The character of the transitions, con- 
tinuous or discontinuous, and their dependence on the 
surface field /3s will be discussed in more detail in a later 
publication. Of course, the transitions between different 
states of a finite chain are no sharp phase transitions, 
as they occur for systems in the thermodynamic limit, 
but rather gradual changes in the weights of various mi- 
crostates of the model system. With increasing chain 
length, the separation between the transition lines de- 
creases so that a phase portrait resembles more closely 
the phase diagram of Fig.[l]as the chain length increases. 
In agreement with theoretical predictions for the ISAWj^ 
we find that the transition line describing the coil-globule 
transition in good solvent is horizontal and perpendicular 
to the adsorption transition line. In contrast to the IS AW 
model, where the crystalline phase always has simple cu- 
bic symmetry and the most compact conformation is cu- 
bic, the BF model supports more than one ordered phase 
and the most compact conformation is highly faceted. 
This leads to a complex sequence of transitions in the 
poor solvent regime, which we will investigate in more 
detail in later work. 



F. Relation to real polymers near surfaces 

Transitions in polymers near surfaces may be induced 
by changing the solvent quality or, more typically, by 
changing the temperature. Since the BF model employed 
in this work orders without forming lamellae and since 
polymer crystallization is typically dominated by kinetic 
effect^2ll a direct comparison with available experimental 
data on surface crystallization of polymers is not possi- 
ble. However, we will discuss some qualitative aspects 
of ordered chains near surfaces. Experiments on alkanes 
physisorbed on graphene surfaces at low coverage show 
the chains to be rod-like (all trans conformations) and 
oriented parallel to the surface.'S^^ Experiments on thin 
films of a variety of crystallizable polymers show that the 
orientation of the chains relative to surface depends on 
the thickness of the layer (see, for example, Ref. [Mjl. 
For very thin films, one typically finds "edge-on" lamel- 
lae corresponding to chains with their backbones ori- 
ented perpendicular to the surface. For example, for 
poly(ethylenc oxide) (PEO) films on bare silicon wafers, 
crystallites with chain backbones perpendicular to the 
surface grow from a layer of adsorbed, non-crystallized 
chains oriented parallel to th e surface through a partial 
dewetting of the surface .'^^'^ 

The fields (3b and /3s employed in this work combine 



the inverse temperature /3 and the interaction parameters 
for net monomer-monomer interactions Sb and monomer 



surface interactions e^, respectively (see Eq. ( 12 1). In our 
phase portrait of Fig. [12] the origin corresponds to the 
limit of infinite temperature while a range of tempera- 
tures at constant interaction parameters corresponds to 
a segment of a straight line through the origin with the 
slope given by the ratio eb/^s of the interaction param- 
eters. For example, for systems without attractive sur- 
face interactions (hard surface) a decrease in temperature 
moves the system up along the vertical line with /?s = 
in Fig. [12] and leads to the two-stage transition from des- 
orbed extended to desorbed compact states. When the 
monomer-monomer interaction parameter et, has a much 
smaller magnitude than the surface interaction parame- 
ter es , a decrease in temperature moves the system along 
an almost horizontal line from (0, 0) to higher /3s values in 
Fig. |12] and leads to the adsorption transition from des- 
orbed (DE) to adsorbed (AE) extended chain conforma- 
tions. For systems where the interaction parameters Cb 
and Es are comparable in size, our phase portrait predicts 
transitions from the disordered DE states to the ordered 
three-dimensional LS states upon lowering the tempera- 
ture. While the details of the LS states clearly depend on 
the model, layered states are not uncommon. Off-lattice 
simulations of confined polymers, for example, show lay- 
ered states for strong confinemenlPSl and the adsorption 
of small molecules from the vapor onto solid substrates 
leads to the occurrence of layered structures ( "multilayer 
adsorption" ) . When the surface interactions are some- 
what larger than the monomer-monomer interactions, a 
decrease in temperature is expected to induce chain ad- 
sorption to the disordered AE states, followed by transi- 
tions to the ordered LS states. During the ordering, the 
chain partially desorbs from the surface to reach the LS 
states. A combination of partial dewetting and ordering 
in a three-dimensional crystallite is also observed in ex- 
periments on PEO on silicon. Finally, when the surface 
interactions are very strong, a transition between two- 
dimensional disordered (AE) and ordered (AC) states is 
expected from our calculations. For real polymeric sys- 
tems, such strong surface attractions might be provided 
by specific interactions (for example hydrogen bonding) 
between the surface and the polymer. 

While some polymers undergo a crystallization transi- 
tion as the temperature is lowered there are many poly- 
mers that remain amorphous. For example, in atactic 
poly(butadiene) the random distribution of monomers of 
different stereoregularity prevents the crystallization of 
the chains. In order to describe non-crystallizable poly- 
mers, one introduces an element of randomness in the 
composition of the model polymer. In this case, we ex- 
pect the phase portrait in field space to be qualitatively 
the same as shown in Fig. [12] for low /3f, values, but to 
contain single lines of continuous transitions between dis- 
ordered extended and compact states. In addition, there 
may be minor transitions in the compact regions corre- 
sponding to local reordering of sections of a chain. 
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V. SUMMARY AND CONCLUSIONS 

In this work we performed Monte Carlo simulations of 
a bond-fluctuation model for a tethered chain with two- 
dimensional Wang-Landau algorithms and umbrella sam- 
pling. The simulations yield density of states results that 
have been evaluated for interaction parameters spanning 
the range from good to poor solvent conditions and from 
repulsive to strongly attractive surfaces. For given fields 
(3s and j3h, we calculated expectation values for chain 
dimensions, contact numbers for surface and monomer- 
monomer contacts, and fluctuations in the contact num- 
bers. Three types of transitions were investigated in some 
detail. For adsorption in good solvent, we compared our 
results with Metropolis Monte Carlo data for the same 
model and found good agreement. For the collapse tran- 
sition in three dimensions, we considered chains tethered 
to a hard surface and found them to behave very simi- 
lar to free chains, with the differences between the two 
situations decreasing with increasing chain length, as ex- 
pected. For the collapse transition in two dimensions, we 
found a dual transition with the same sequence of tran- 
sitions that is observed for three-dimensional chains: a 
coil-globule transition that changes the overall chain size 
is followed by a local rearrangement of chain segments. 

In order to investigate the overall phase behavior of the 
tethered chains considered in this work, we located max- 
ima of the susceptibilities Xs and Xb in the Ps-Pb plane. 
We found that all transitions from good to poor solvent 
conditions are dual in nature and that the separation be- 
tween the two lines belonging to the same transition de- 
creases with increasing chain length. In agreement with 
theoretical predictions for the IS AW model, * we find that 
the transition line describing the coil-globule transition 
in good solvent is horizontal and perpendicular to the 
adsorption transition line. For poor solvent conditions, 
early work on the ISAW model predicted the exi stence 
of a surface attached globule (SAG) phas^^^^^^^ which 
was later found to consist of a whole sequence of layered 
states.!^ In the corresponding parameter region we find 
one transition (between single and double layer states) 
that is similar to the AC to LS (SAG) transition of the 
ISAW. The next transition of the BF model, however, 
changes the symmetry of the ordered segments and has 
no correspondence in the ISAW model. From the max- 
ima of the susceptibilities of finite-length chains alone it 
is not possible to resolve the nature of the transitions 
near points where transition lines meet. This is an in- 
teresting question that we will try to address in future 
work. 
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APPENDIX A: SIMULATION DETAILS 

In this section, we present details of the density-of- 
states and production simulations for chains of length 
iV = 16, TV = 32, and N = 64. For each chain 
length, density of states results were generated with 
Wang-Landau algorithms and then refined with umbrella 
sampling. In order to compare results from groups of sim- 
ulations, the average over the ln{g) values in the interior 
(7V/4 < < 37V/4, 57V/16 < Ub < 257V/16) of the range 
of accessible states is calculated for each simulation and 
the lii(g) values are shifted by this amount. To deter- 
mine uncertainties for the WL results, the average and 
standard deviation of the ln{g) values are calculated for 
each state {ns,nb). For umbrella sampling results, the 
log-density of states values are weighted by the length of 
the umbrella sampling simulation and the weighted av- 
erage ln{g) and standard deviation crin(g) are calculated. 
Parameters for the umbrella sampling simulations and 
some numerical characteristics of the final densities of 
states are presented in Table |T] In the following, we high- 
light differences in the simulations for the different chain 
lengths. Production was carried out either in separate 
simulations or concurrently with the construction of the 
density of states. A summary of the production simula- 
tion parameters is presented in Table [ll] 

TABLE I: Umbrella sampling parameters and some density 
of states characteristics for chains of length A'^ = 16, A'^ = 32, 
and — 64. A'^„ is the number of results for which um- 
brella sampling simulations were performed, Ku is the length 
of the simulations in Monte Carlo steps, Ns is the number of 
states sampled, ni,,max(l) and ni,,niax(A'') represent the largest 
number of bead-bead contacts for Us = 1 and Us = N sur- 
face contacts, respectively. The values for the range of the 
log-density of states {A\n{g)), and its uncertainties {(Jin(g)) 
represent results after umbrella sampling. For A'^ = 16 and 
N — 32, the states in this table are believed to be the com- 
plete set; for A' — 64, only states included in our evaluation 
are represented. 



TV 


16 


32 


64 




4 


3 


2 




1 


16, 10, 18 


18.5, 18.8 




748 


3467 


15268 


^6, max ( 1) 


53 


128 


288 


Jib, max (Af) 


34 


76 


163 


Aln(g) 


50.8 


109.5 


227 


average oi^f^g^ 


0.008 


0.008 


0.015 


median o"in(a) 


0.004 


0.003 


0.004 


maximum ain(t,) 


0.18 


0.5 


1.5000 



TV = 16 For chains of length TV = 16, four simula- 
tions with the original Wang-Landau algorithm were per- 
formed. The histogram was considered flat when the vis- 
its to each individual (n^, Uh) state were no less than 0.8ft., 
where h is the average of the visits to all states. The sim- 
ulations required between 1.2x10^ and 1.9x10® Monte 
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TABLE II: Production parameters for chains of length = 
16, A'^ — 32, and A = 64. The table entries represent the 
simulation length in MC steps, where configurations are eval- 
uated every 10 MC steps. The left column indicates the type 
of simulation. For simulations sampling with the density of 
states, the number of replicas (Nr) is indicated and, in one 
case, the maximum number of bead contacts considered. For 
simulations with the Metropolis acceptance criterion, the val- 
ues for the field variables Ps and /3b are shown. 



N 


16 


32 


64 


Nr 


= 1 




1 X 10** 


4 X 10^ 


3.35 X 10'' 


Nr 


= 3 




1 X 10^ 






Nr 


= 4, 


nt < 256 






2.6 X 10** 


Nr 


= 5 






1 X 10^ 


3 X 10** 


Ps 


= /3i> 


= 


1 X 10** 


1 X 10^ 


1.5 X 10^ 


f3s 


^Pb 


= -1 






1 X 10^ 


Ps 


= Pb 


= -2 


2 X 10^ 


1 X lO'' 


4 X 10** 



Carlo steps to converge. The results agreed well with 
each other; the average, median, and maximum standard 
deviation of the ln{g) values are 0.017, 0.008, and 0.63, 
respectively. The largest deviations occur for states with 
the largest number of bead contacts for given number of 
surface contacts, i.e. for states near the upper rim of the 
range of accessible states in Fig. |2] Umbrella sampling 
simulations with 10® Monte Carlo steps were performed 
for each of the WL density of states results for = 16 
and evaluated as described in Sec. |IID| The agreement 
between the four different results is excellent, as may be 
seen from the values of the standard deviations presented 
in Table |l] As before, the largest deviations occur at the 
upper rim of the range of accessible states. 

N — 32 For chains of length A^ = 32, simulations with 
the original Wang-Landau algorithm would not converge. 
We modified the flatness criterion and considered aver- 
ages over groups of states. Each group consisted of A^ 
states and two sets of groups were created. One set was 
formed by going up along columns in the {ns,nt,) plane 
and adding states to a group until it was filled, continu- 
ing in the next column when the top was reached before 
the group was complete. The other set was formed by 
going along rows in the same way. The histogram was 
considered flat when the visits to each group of states 
in the two sets was no less than 0.8 and no more than 
1.25 of the average of the visits over all groups. We ob- 
tained results from two simulations (A and B) with a 
Wang-Landau algorithm where the density of states was 
updated only for accepted moves (see Eq. (|3|) and from 
one simulation (C) with a five-replica Wang-Landau al- 
gorithm where the density of states was updated after 
accepted as well as after rejected moves, Eqs. ^ and 
Q. The results for the density of states from the five 
replicas agreed well in the regions of overlap and were 
combined with a simple linear switching function. 

The N = 32 simulations took between 1.3 x 10^° and 



4.1 X 10"'^° MC steps to converge. All simulations missed 
or severely underestimated the density of states values of 
four states on the upper rim of accessible states in Fig. |2] 
Except for states on the rim, the two single-replica results 
A and B agree reasonably well with each other, the av- 
erage, median, and maximum standard deviation of the 
ln(g) values are 0.029, 0.011, and 4.4, respectively. How- 
ever, there are systematic deviations between the l'a{g) 
values for high n(,-states between result C from the five- 
replica simulation and results A and B. 

In order to prepare for umbrella sampling, the density 
of state results for A^ = 32 from the Wang-Landau simu- 
lations were augmented by assigning reasonable guesses 
to the la(g) values that were missed or nearly missed in 
the simulations. We performed umbrella sampling simu- 
lations with the parameters given in Table [l] and updated 
the density of states results as described in Sec. IIP The 
systematic deviations between results A, B and C disap- 
peared after umbrella sampling and the density of states 
values agree very well, (see the values for the standard 
deviations in Table |l]) except for states very close to the 
upper rim of the range of accessible states in Fig. |2] For 
results A and B, umbrella sampling led to significant in- 
creases in the normalized lii(g) values for a large range of 
states with high numbers of bead contacts. For result C, 
umbrella sampling decreased the normalized ln(g) values 
slightly for states with very high numbers of bead con- 
tacts and changed the values significantly on the rim of 
accessible states. One state was not visited in the C sim- 
ulation and its density of state value was calculated from 
the A and B results only. Production simulations were 
carried out with the parameters in Table |llj For some 
states near the rim, production data were supplemented 
by evaluating stored chain conformations. 

TV = 64 For chains of length A^ = 64, one density of 
states result with a cutoff in the number of bead con- 
tacts at nil — 4N was obtained with the modified Wang- 
Landau algorithm described for chains of length A^ = 32 
(simulations A and B), above. This simulation was com- 
pleted in 2 X 10® MC steps; simulations without cutoff 
would not converge. 

In order to obtain density of states results for a larger 
range of bead contact values, we applied the global up- 
date algorithm described in Section [IIC| For the results 
presented here, we used the following shift parameters in 
Eq. (§, K = A = = 10Mn(/fc), where ln(/fc) = 2-('=-i) 
is the increment of the log density of states at refinement 
level k. The first global update simulation used a value 
of p = 10^1n(/fe) for the uniform growth criterion. It 
was started at the fc — 15 level and allowed to progress 
through about 8 x 10^ MC steps before the refinement 
level was set to fc = 18. At this level, about 7.5 x 10* 
MC steps were carried out before our groups-of-state flat- 
ness criterion was satisfied and the simulation progressed 
through the final two levels in a standard Wang-Landau 
algorithm with our fiatness criterion, which took about 
3.5 X 10* MC steps, for a total simulation time of about 
1.9 X 10® MC steps. The density of states values obtained 
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in this global update simulation agreed reasonably well in 
the range of common states with the results obtained in 
the simulation with cutoff described above. It extended 
the range of visited states to higher bead contacts but 
it missed a number of states with large surface and bead 
contact values that had been found in the simulation with 
cutoff. After removing states whose density of states val- 
ues had been severely underestimated, we combined the 
results from the two simulations with a linear switching 
function; this combined density of states is referred to as 
result A in the following. 

The second global update simulation was started at 
the fc = 18 level and used a value of p = 101n(/fe). With 
these parameters, the density of states in the interior 
is generated through a larger number of smaller shifts. 
The simulation was allowed to proceed through a total 
of about 3.15 x 10^*^ MC steps, after which time we col- 
lected the density of states, our result B. During this 
time, global updates and the addition of new states be- 
came increasingly rare. For example, only 31 of the more 
than 15000 visited states were found in the second half 
of the simulation. The range of visited states from this 
simulation is larger than that of result A. The differences 
between results A and B are largest for states with very 
high bead-contact values; the average, median, and max- 
imum standard deviation of the ln(g) values are 0.16, 0.1, 
and 18.3. In order to prepare for umbrella sampling, the 
density of state results A and B were smoothed by as- 
signing reasonable guesses to the ln{g) values that were 
missed or nearly missed in the simulations. 

Unfortunately, the umbrella sampling simulations for 



both A and B missed some states completely (20 for A 
and 12 for B out of a total of 15275 input states), nearly 
missed a few others (8 for A and 2 for B), and found 
some states for the very first time (12 for A and 14 for 
B, with 4 in common for a total of 22 new states). All 
of the problematic states are very close to the upper rim 
of the range of accessible states shown in Fig. |2] Except 
for the newly found states, which were discarded, the 
In (5) values were updated at the end of the umbrella 
sampling as described in Section |IID[ For states that 
were found in both umbrella sampling simulations, the 
log-density of states results were combined as described 
above. For states that had been missed or nearly missed 
in one of the umbrella sampling simulations, we used the 
density of states result from the other simulation and 
assigned an uncertainty of cin(g) = 1-5. Except near the 
rim, the agreement between the results is very good (see 
Table |l]) . Production simulations were carried out with 
the parameters in Table |TT] For some states near the rim, 
production data were supplemented by evaluating stored 
chain conformations. For 13 states on the rim, values of 
the quantities of interest were estimated by extrapolation 
from data with smaller numbers of bead contacts and 
generous errors were assigned. The impact of the rim 
states depends on the values of the field variables, in 
particular on the value of /?;,, as discussed in Section [lV A[ 
In this work, we focus on conditions where their impact is 
small. For future evaluation, we are currently performing 
simulations with a slightly modified simulation algorithm 
that we hope will improve the density of states results 
near the rim. 
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